{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Angle-distribution function\n",
    "\n",
    "The angle-distribution function describes the probability density of finding a triple of atoms with a specific angle. Note that this requires restricting the distance between atoms to a specific bond-length, or some other means to find bonds. Here we describe how to compute the angle distribution if bonds can be determined from a simple distance-based criterion, i.e. two atoms are considered bonded if their distance is smaller than $r_\\textrm{b}$. The angle distribution is then\n",
    "\n",
    "$$\n",
    "f_\\text{a}(\\phi) = \\frac{1}{N} \\sum_{i=1}^N \\frac{1}{N_i^2} \\sum_{j=1}^N \\sum_{k=1}^N \\theta(r_{ij}-r_\\textrm{b}) \\theta(r_{ik}-r_\\textrm{b})\\, \\delta\\left(\\text{arccos}(\\vec{r}_{ij}\\cdot \\vec{r}_{ik}/(r_{ij} r_{ik}))-\\phi\\right)\n",
    "$$\n",
    "\n",
    "where $\\theta(x)$ is the Heaviside step function, $N$ is the total number of atoms and $N_i=\\sum_j \\theta(r_{ij}-r_\\textrm{b})$ is the number of neighbors (coordination number) of atom $i$. This angle-distribution function has the property $\\int \\text{d}\\phi\\, f_\\text{a}(\\phi)=1$, i.e. it is a probability density. It can be calculated from using the utility `angle_distribution` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 3.141592653589793)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAG0CAYAAAAy8S2PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8WgzjOAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYu0lEQVR4nO3deVyU9fo//tfMwDCA7DuKgEseDQUjUdDSjuaWZvWt0/GUmmWf7Gelknm008dK61CW2GZy1NxOdVrc+pwWE8mlFFckUXMF2WTfd4aZ+/fHMAOjgDD3wGyv5+MxjyP33HPfF/dBvLre7/f1lgiCIICIiIjIBkhNHQARERFRT2HiQ0RERDaDiQ8RERHZDCY+REREZDOY+BAREZHNYOJDRERENoOJDxEREdkMO1MHYG7UajVu3LgBFxcXSCQSU4dDREREnSAIAqqqqhAYGAiptP26DhOfm9y4cQNBQUGmDoOIiIgMkJ2djT59+rT7PhOfm7i4uADQPDhXV1cTR0NERESdUVlZiaCgIN2/4+1h4nMT7fCWq6srEx8iIiILc7tpKpzcTERERDaDiQ8RERHZDCY+REREZDOY+BAREZHNYOJDRERENoOJDxEREdkMJj5ERERkM5j4EBERkc1g4kNEREQ2g4kPERER2QwmPkRERGQzmPgQERGRzWDiQ0RkBeqVKlOHQGQRmPgQEVm4LUcy8Kf/3YvECwWmDoXI7DHxISKyYCq1gA2H0wEA25OvmzYYIgvAxIeIyIIdvVaMvIr65j+XoKS6wcQREZk3Jj5ERBZsx+kc3Z9VagE/n+dwF1FHmPgQEVmoynol9p7LBwBMCfMHAPyYlmfKkIjMHhMfIiIL9cPZPDQ0qTHQtxeWTxkMQDP0xeEuovYx8SEislDaYa5HI/ugr5cThvZ2g1oAh7uIOsDEh4jIAl0rqsbpzDJIJcDDw3sDAB4YFgAA+CHthilDIzJrTHyIiCzQzuZqz9g7fODrqgAAPDBUk/gkXytBMYe7iNrExIeIyMKo1AJ2peQCAB6NDNIdD/J0wrA+2uGufFOFR2TWmPgQEVmYI1eLkV9ZDzdHe4wf7Kv33tTmqs8PZ7m6i6gtTHyIiCyMdlLzg+GBUNjL9N7TDncdS+dwF1FbmPgQEVmQijqlbhjr0cg+t7wf5OmE8ObhLm2PHyJqwcSHiMiCtO7dM6yPW5vncLiLqH1MfIiILMiO09kANNUeiUTS5jnaxOd4RgmKqjjcRdQaEx8iIgtxragaKVnlkEklut49bdEb7uLqLiI9THyIiCxEW7172qNtZvgjh7uI9DDxISKyAPq9e26d1Hyz1sNdhVX13RobkSVh4kNEZAE66t3Tlj4eTggPctc0M+TqLiIdJj5ERBZA27tnRkQgHOxktzlbY5p2dVcah7uItJj4EBGZudv17mnPlKH+AIDjGaUc7iJqZraJT1xcHEaMGAEXFxf4+vrioYcewqVLlzr8zNatWyGRSPReCkXHEwCJiMydtnfPHX69MLR327172tLHwwkRQe4Q2MyQSMdsE59Dhw5hwYIFOHbsGBITE6FUKjFx4kTU1NR0+DlXV1fk5eXpXpmZmT0UMRFR9+hM7572TBvGZoZErdmZOoD27N27V+/rrVu3wtfXF6dPn8a9997b7uckEgn8/f27Ozwioh7RunfPQxHt9+5pz5ShAXjrhz9w4nopCivrb7sMnsjamW3F52YVFRUAAE9Pzw7Pq66uRnBwMIKCgjBjxgycP3++w/MbGhpQWVmp9yIiMhdd6d3Tlt7ujhjet3m4i80MiSwj8VGr1Vi0aBFGjx6NsLCwds8bNGgQNm/ejO+++w6ff/451Go1YmJikJOT0+5n4uLi4ObmpnsFBQV1x7dARNRlXe3d0x7tju3fc7iLCBJBEARTB3E7zz//PH766Sf89ttv6NOn83/5lUolBg8ejJkzZ2LVqlVtntPQ0ICGhpa9bCorKxEUFISKigq4urqKjp2IyFCHLxdh9uYTcHO0x4l/jO/0Mvab3SivQ8w7v0AiAY4vH8/hLrJKlZWVcHNzu+2/32Zf8XnhhRfw/fff48CBA11KegDA3t4ew4cPx9WrV9s9x8HBAa6urnovIiJzcOhyEQBNF2ZDkx4ACHR3xF3Nw10/cXUX2TizTXwEQcALL7yA3bt345dffkFoaGiXr6FSqZCWloaAgIBuiJCIqHtdL9asYh0SKP4/yLRbWHB1F9k6s018FixYgM8//xxffvklXFxckJ+fj/z8fNTV1enOmT17NpYvX677euXKldi3bx/S09ORkpKCJ598EpmZmZg3b54pvgUiIlEyS2sBACFeTqKvpU18TmaWoqCSzQzJdplt4rN+/XpUVFRg3LhxCAgI0L2+/vpr3TlZWVnIy2v5r5eysjI8++yzGDx4MKZOnYrKykocPXoUQ4YMMcW3QERkMJVaQFaJNvFxFn291sNdBy8Vir4ekaUy2z4+nZlzffDgQb2v165di7Vr13ZTREREPSe/sh6NKjXspBIEuBlnMnJ4kDtSsspxtbDaKNcjskRmW/EhIrJlmc3ze4I8nWAnM86v6n4+vQAAGcUdd8AnsmZMfIiIzJB2fk+wEeb3aPXz1gyZpRcx8SHbxcSHiMgMXS/RJCfGmN+jFdqc+GSV1kKpUhvtukSWhIkPEZEZyizWVHz6ehqv4uPvqoCjvQxNagE5ZXW3/wCRFWLiQ0RkhnQVH2/jJT5SqQQhuuEuTnAm28TEh4jIzAiCgCzdHB/jDXUBLfN8OMGZbBUTHyIiM1NU3YDaRhWkEqCPh6NRr93Pp7niw8SHbBQTHyIiM5PZ3Lgw0N1R1B5dbQnlUBfZOCY+RERmRrtHlzGXsmuFcqiLbBwTHyIiM6Ot+Bh7fg8A9PPWNDEsqGxATUOT0a9PZO6Y+BARmRljbk56Mzcne3g5ywGw6kO2iYkPEZGZySzRDnUZv+IDtJrnw8SHbBATHyIiMyIIgq4SY8yuza1pV3ZlcOsKskFMfIiIzEh5rRJV9Zq5N8bs2txaaPM8n/Riruwi28PEh4jIjGg7Nvu5OsBRbtyl7Fpc2UW2jIkPEZEZ6a6Oza31bzXUJQhCt92HyBwZJfFRKpXIzs7GpUuXUFpaaoxLEhHZpOvF3beiS6uvlxMkEqCqoQlF1Q3ddh8ic2Rw4lNVVYX169dj7NixcHV1RUhICAYPHgwfHx8EBwfj2WefxcmTJ40ZKxGR1evuFV0A4GAn022FwQnOZGsMSnzi4+MREhKCLVu2YMKECdizZw9SU1Nx+fJlJCcn4/XXX0dTUxMmTpyIyZMn48qVK8aOm4jIKl0v6b6uza1pGxlyng/ZGjtDPnTy5EkcPnwYd955Z5vvR0VF4emnn0ZCQgK2bNmCX3/9FQMHDhQVKBGRLcjSNS/svooPoJngfOhyEXv5kM0xKPH5z3/+06nzHBwcMH/+fENuQURkc6rqlSiubgSgmYfTnXS7tHOoi2yM0VZ1HT9+3FiXIiKySdo9uryc5XBV2HfrvVqGutjLh2yL0RKfxx57zFiXIiKySdrEp7urPQAQ2lzxySqtRZNK3e33IzIXXRrq+stf/tLmcUEQuIydiEgk7cTm7p7fAwABrgoo7KWoV6qRU1aHEO/uvyeROehS4rN//378+9//Rq9evfSOC4KAw4cPGzUwIiJbk1WibV7Y/RUfqVSCEC9nXMyvQnpxNRMfshldSnzGjRsHFxcX3Hvvvbe8N2zYMKMFRURki3qy4gNoJjhfzK9CelEN/vynHrklkcl1KfHZtWtXu+8lJiaKDoaIyJZl9mDFB+CeXWSbRE1uzs/PN1YcREQ2ra5RhfzKegDd27W5Ne3KLi5pJ1siKvGZOHGiseIgIrJp2saFLgo7eDh171J2Le3KLlZ8yJaISny4qy8RkXFktprfI5FIeuSe/ZqHuvIr61HT0NQj9yQyNVGJT0/95SQisnY9Pb8HANyd5PB0lgNg1Ydsh9EaGBIRkeF6ekWXFic4k61h4kNEZAZ6smtza/2Y+JCNEZX4yGQyY8VBRGTTMktNVPHRbVbKPbvINohKfM6cOWOsOIiIbFZjkxq5ZXUAgBBWfIi6FYe6iIhMLKesFmoBcLSXwcfFoUfv3c+nuZdPcQ1X6pJNEJ341NXVoba2Vvd1ZmYmPvjgA+zbt0/spYmIbELrFV09vVq2r6cTJBKgqr4JxdWNPXpvIlMQnfjMmDED27dvBwCUl5dj5MiRWLNmDWbMmIH169eLDpCIyNppV3T15FJ2LYW9DL3dHQFwuItsg+jEJyUlBffccw8AYMeOHfDz80NmZia2b9+Ojz76SHSARETWTlvx6emJzVra4a6MYk5wJusnOvGpra2Fi4sLAGDfvn145JFHIJVKMWrUKGRmZooOkIjI2mXqKj4mSny8tSu7WPEh6yc68RkwYAD27NmD7Oxs/Pzzz7r9uwoLC+Hq6io6QCIia9dS8en5oS6gpYlhOoe6yAaITnxWrFiBJUuWICQkBCNHjkR0dDQATfVn+PDhogMkIrJmKrWA7DLTNC/U6sfNSsmG2Im9wKOPPooxY8YgLy8P4eHhuuPjx4/Hww8/LPbyRERW7UZ5HZQqAXKZFAFujiaJQVvxySypQZNKDTsZO52Q9RKd+ACAv78//P399Y5FRUUZ49JERFZNO8wV5OkImdQ0Gz8HujnCwU6KhiY1csvrTDbXiKgnGCXxSUpKQlJSEgoLC6FWq/Xe27x5szFuQURklUy1OWlrUqkEod7OuJhfhfTiGiY+ZNVE1zPffPNNTJw4EUlJSSguLkZZWZnei4iI2mfqFV1aoVzZRTZCdMUnISEBW7duxaxZs4wRDxGRTbneqmuzKYXq9uxiLx+ybqIrPo2NjYiJiTFGLERENifLTBKfliaGrPiQdROd+MybNw9ffvmlMWIhIrIparWAzFLTz/EBONRFtkP0UFd9fT02bNiA/fv3Y9iwYbC3t9d7Pz4+XuwtiIisUmFVA+qVasikEvT2MM1Sdi1t9+a8inrUNjbBSW6UtS9EZkf0T/bZs2cREREBADh37pzeez29yzARkSXRrujq7e4IexP3zvFwlsPDyR5ltUpkFNfgzkA3k8ZD1F1EJz4HDhwwRhy3iIuLw65du3Dx4kU4OjoiJiYG7777LgYNGtTh57799lv87//+L65fv46BAwfi3XffxdSpU7slRiIiMcxlfo9WqLczyrLKmfiQVTNaLfPChQvIyspCY2Oj7phEIsH06dMNut6hQ4ewYMECjBgxAk1NTXj11VcxceJEXLhwAc7ObY+FHz16FDNnzkRcXBymTZuGL7/8Eg899BBSUlIQFhZmUBxERN3FHHr4tBbq3QspWeXI4DwfsmKiE5/09HQ8/PDDSEtLg0QigSAIAFqGuVQqlUHX3bt3r97XW7duha+vL06fPo177723zc98+OGHmDx5Ml555RUAwKpVq5CYmIhPPvkECQkJBsVBRNRdMs2s4qPds4ublZI1Ez2ovHDhQoSGhqKwsBBOTk44f/48Dh8+jLvvvhsHDx40QogaFRUVAABPT892z0lOTsaECRP0jk2aNAnJycntfqahoQGVlZV6LyKinmBuFZ9+3KWdbIDoxCc5ORkrV66Et7c3pFIppFIpxowZg7i4OLz00kvGiBFqtRqLFi3C6NGjOxyyys/Ph5+fn94xPz8/5Ofnt/uZuLg4uLm56V5BQUFGiZmIqCOCIJhdxSdUu0t7UbWuek9kbUQnPiqVCi4uLgAAb29v3LhxAwAQHByMS5cuib08AGDBggU4d+4cvvrqK6Ncr7Xly5ejoqJC98rOzjb6PYiIblZa04jqhiZIJECQp3kkPiFezpBIgMr6JpTUNN7+A0QWSPQcn7CwMPz+++8IDQ3FyJEjsXr1asjlcmzYsAH9+vUTHeALL7yA77//HocPH0afPn06PNff3x8FBQV6xwoKCm7ZOb41BwcHODg4iI6TiKgrtFtVBLgqoLCXmTgaDYW9DIFujsgtr0NGcQ28e/F3I1kf0RWf1157Tbcj+8qVK5GRkYF77rkHP/74Iz766CODrysIAl544QXs3r0bv/zyC0JDQ2/7mejoaCQlJekdS0xMRHR0tMFxEBF1B3PZnPRm/XTDXZznQ9ZJdMVn0qRJuj8PGDAAFy9eRGlpKTw8PEQ1MFywYAG+/PJLfPfdd3BxcdHN03Fzc4Ojo6bD6ezZs9G7d2/ExcUB0Ey0Hjt2LNasWYMHHngAX331FU6dOoUNGzaI+A6JiIxPW/EJ8TaPYS6tft7O+PVKMa5xs1KyUt3SKtTT01N01+b169ejoqIC48aNQ0BAgO719ddf687JyspCXl6e7uuYmBh8+eWX2LBhA8LDw7Fjxw7s2bOHPXyIyOxkNVd8+nqaV8VHt0s7Kz5kpQyq+MTGxmLVqlVwdnZGbGxsh+cauldXZ1YUtLVc/rHHHsNjjz1m0D2JiHqKruJjJiu6tEK5SztZOYMSnzNnzkCpVOr+3B7u1UVE1DaznePTXPHJLKmFSi1AJuXvcbIuBiU+rffn6q69uoiIrFVFrRJltZr/eDSXHj5age6OkNtJ0dikRm5ZHfqaWXxEYpl2O2AiIhuUWaqp9nj3coCzg9G2TDQKmVSC0OYqFCc4kzUyeI5PZxk6x4eIyFplmun8Hq0QbydcKqhCZnENMMjU0RAZl8FzfFpLSUlBU1MTBg3S/A25fPkyZDIZIiMjxUdIRGRlsss0iU9fM+nYfDPt3mHaCdhE1kT0HJ/4+Hi4uLhg27Zt8PDwAACUlZVh7ty5uOeee4wTJRGRFSmsbAAA+LkpTBxJ27QTrrUTsImsieg5PmvWrEFcXJwu6QEADw8PvPXWW1izZo3YyxMRWZ2CynoAgJ+LeW4JoR2Cy2TFh6yQ6MSnsrISRUVFtxwvKipCVVWV2MsTEVmdwqrmio+rmVZ8mpe0Z5fVokmlNnE0RMYlOvF5+OGHMXfuXOzatQs5OTnIycnBzp078cwzz+CRRx4xRoxERFZFW/HxdTXPik+AqwJyOymUKgF5FfWmDofIqESvo0xISMCSJUvwt7/9DUqlEoIgwN7eHs888wzee+89Y8RIRGQ1BEHQzfHxdTHPio9UKkFfTydcLazG9ZIaBJnpJGwiQ4iu+Dg5OeHTTz9FSUkJzpw5g9TUVJSWluLTTz+Fs7N5dSQlIjK1ijolGpuHj8y14gO0zPPhyi6yNkbrnJWZmYkbN26gsbER169f1x1/8MEHjXULIiKLV9Bc7XF3soeDnczE0bRPt7KLe3aRlRGd+KSnp+Phhx9GWloaJBKJbnNR7T5dKpVK7C2IiKxGy4ou8xzm0mLFh6yV6KGuhQsXIjQ0FIWFhXBycsK5c+dw+PBh3H333W3unk5EZMu0K7rMeZgLYC8fsl6iKz7Jycn45Zdf4O3tDalUCplMhjFjxiAuLg4vvfRSh7u3ExHZGt2KLrOv+DQnPqW1UKsFSLlLO1kJ0RUflUoFFxcXAIC3tzdu3LgBAAgODsalS5fEXp6IyKoUaoe6zLziE+iugJ1UgsYmNfIruaSdrIfoik9YWBh+//13hIaGYuTIkVi9ejXkcjk2bNiAfv36GSNGIiKrYe7NC7XsZFIEeToho7gG10tqEOjuaOqQiIxCdMXntddeg1qtWZq5cuVKZGRk4J577sGPP/6Ijz76SHSARETWpGWoy7wrPgAQzK0ryAqJrvhMmjRJ9+cBAwbg4sWLKC0thYeHh25lFxERaWiXs/uaecUH0M7zKcJ1TnAmKyKq4qNUKjF+/HhcuXJF77inpyeTHiKimwiCgCLdUJcFVXyKWfEh6yEq8bG3t8fZs2eNFQsRkVUrr23p2uxjAUNd2pVdrPiQNRE9x+fJJ5/EZ599ZoxYiIisWkGVZn6Ph5l3bdZqPcdH25yWyNKJnuPT1NSEzZs3Y//+/YiMjLxlf674+HixtyAisgrazUnNfUWXVh8PJ0glQJ1ShaKqBouYl0R0O6ITn3PnzuGuu+4CAFy+fFnvPc7zISJqoV3RZQnDXAAgt5Oit4cjskvrcL2klokPWQXRic+BAweMEQcRkdWzlB4+rYV4OTcnPjWICvU0dThEohk0xycrK6tL5+fm5hpyGyIiq2IpXZtba5nnwwnOZB0MSnxGjBiB5557DidPnmz3nIqKCmzcuBFhYWHYuXOnwQESEVkLXQ8fM9+nq7WWlV1c0k7WwaChrgsXLuDtt9/G/fffD4VCgcjISAQGBkKhUKCsrAwXLlzA+fPncdddd2H16tWYOnWqseMmIrI42lVdllXx4S7tZF0Mqvh4eXkhPj4eeXl5+OSTTzBw4EAUFxfrGhk+8cQTOH36NJKTk5n0EBE1K7Sgrs1aIa2aGHJJO1kDUZObHR0d8eijj+LRRx81VjxERFZJEAQUVlnOPl1aQZ5OkEiAqoYmlNY0wquX5cRO1BbRDQyJiOj2ymqVUKo0FRNLWc4OAAp7GQKaK1Sc50PWgIkPEVEP0FZ7PJ3lFtG1uTXO8yFrwsSHiKgHtKzospxqj1aIt2aeDys+ZA2Y+BAR9QBt12ZLmtisxYoPWRMmPkREPaBI27XZEis+Xqz4kPUQvWUFACiVSuTn56O2thY+Pj7w9GRbcyKi1loqPpaX+LDiQ9bE4IpPVVUV1q9fj7Fjx8LV1RUhISEYPHgwfHx8EBwcjGeffbbDzs5ERLakQLddhSUOdWkqPuW1SpTXNpo4GiJxDEp84uPjERISgi1btmDChAnYs2cPUlNTcfnyZSQnJ+P1119HU1MTJk6ciMmTJ+saGxIR2SrtBqWWtF2FlpPcTjcpm8NdZOkMGuo6efIkDh8+jDvvvLPN96OiovD0008jISEBW7Zswa+//oqBAweKCpSIyJK1dG22vKEuQLNnV2FVAzJLahAR5G7qcIgMZlDi85///KdT5zk4OGD+/PmG3IKIyGq07tpsiUNdgGa468T1UlwvZsWHLJtRJjefPXsWv/76K+RyOWJiYtqtBBER2SK9rs0WuuVDiDcnOJN1EJ34fPjhh1i8eDFcXV0hk8lQVlaGoUOHYtu2bYiIiDBCiERElk07sdnTWQ65nWV2EQlpXtl1nYkPWTiD/gZu3rwZKSkpaGhowNtvv4133nkHZWVlKCkpQXp6OqZMmYJ77rkHR48eNXa8REQWR7eU3QJ7+GhpV3ZlcnIzWTiDKj7vv/++bqWWWq3GyZMn8eGHH2L48OGIiIjAO++8g6CgICxZsoTJDxHZPO3EZkud3wO0JD4lNY2orFfCVWFv4oiIDGNQxefChQuoqqrC0aNHYW9vD6lUiq+++gpTp06Fp6cn+vXrh927d+P06dP44YcfcP36dSOHTURkObQTmy254uOisId3LzkAIItVH7JgBg82KxQKjBgxAqNHj0Z4eDiOHTuGqqoqpKWl4a233sKAAQOgVCoxe/Zs9OvXD66ursaMm4jIYhRYQcUHaOngzHk+ZMlET25es2YNxo0bh/T0dMyfPx/h4eEICgpCSkoKAgMDkZOTg5ycHJw7d84Y8RIRWZyWrs2WW/EBNMNdpzPLOM+HLJroxCciIgKnT5/G/PnzMWrUKAiCZsmmnZ0dNm/eDADo06cP+vTpI/ZWREQWSdu12ccCuza3plvZVcyKD1kuo/Tx6d+/PxITE1FQUIBjx46hsbER0dHRTHaIiAAUWlHFB+DKLrJsRkl8tPz8/DBjxgxjXpKIyKKp1YKu4mPpc3zYy4esgUGTm7Oysrp0fm5uriG3ISKyeGW1jWhSa6YAeFto12YtbeJTWNWA2sYmE0dDZBiDEp8RI0bgueeew8mTJ9s9p6KiAhs3bkRYWBh27txpcIBERJZMu6LLy4K7Nmu5OdnD3UnTv4fDXWSpDO7j4+zsjPvvvx/+/v544IEH8Oyzz+LFF1/Ek08+ibvuugu+vr7YvHkzVq9ejZdeesmg4A4fPozp06cjMDAQEokEe/bs6fD8gwcPQiKR3PLKz8836P5ERGIVaHv4WPgwl5Z2STv37CJLZVDi4+Xlhfj4eOTl5eGTTz7BwIEDUVxcrOvm/MQTT+D06dNITk7G1KlTDQ6upqYG4eHhWLduXZc+d+nSJeTl5elevr6+BsdARCRGUXPFx5KbF7YW0jzB+TorPmShRE1udnR0xKOPPopHH33UWPHomTJlCqZMmdLlz/n6+sLd3d34ARERdZG19PDRYsWHLJ1RV3W1plQqkZWVhYaGBt2xIUOGdNft9ERERKChoQFhYWF44403MHr06HbPbWho0IuxsrKyJ0IkIhuhHeqy9BVdWrqKTzErPmSZumWm3X/+8x9ERERg+PDhmDNnDiIiIjBv3rzuuJWegIAAJCQkYOfOndi5cyeCgoIwbtw4pKSktPuZuLg4uLm56V5BQUHdHicR2Y5CKxvqYsWHLF23JD5xcXE4efIk+vXrh5MnT+LEiRM9klAMGjQIzz33HCIjIxETE4PNmzcjJiYGa9eubfczy5cvR0VFhe6VnZ3d7XESke0oaO7hYy2Tm7UVnxsV9ahXqkwcDVHXdUvi4+DgACcnzV8OpVKJiIgInD9/vjtudVtRUVG4evVqu+87ODjA1dVV70VEZCwtXZutI/HxdJbDxUEzSyK7lMNdZHkMSnyWLl2K+vr6dt/39/dHeXk5pk+fjilTpuDxxx832RBSamoqAgICTHJvIrJtarWAoirrGuqSSCQI9ubKLrJcBk1u/uCDD/Dyyy9DoVDgqaeewqeffqqr8ADAf//7XwDAqlWrcPDgQVRWVmLy5Mldvk91dbVetSYjIwOpqanw9PRE3759sXz5cuTm5mL79u26uEJDQ3HnnXeivr4emzZtwi+//IJ9+/YZ8m0SEYlS2qprs4+VJD6AZp7PudxKzvMhi2RQxScwMBCpqakAgH//+9+orq5u99xx48bhwQcfhFwu7/J9Tp06heHDh2P48OEAgNjYWAwfPhwrVqwAAOTl5eltn9HY2IiXX34ZQ4cOxdixY/H7779j//79GD9+fJfvTUQklnYpu3cvOexllt21ubWWXj5MfMjyGFTxefnllzF9+nSMHDkSAPDFF19g9OjRGDp0KBwdHY0W3Lhx4yAIQrvvb926Ve/rpUuXYunSpUa7PxGRGNrNSX1crGN+j1bLyi4OdZHlMeg/QV588UWcOnUKkydPhiAIWLduHWJiYuDq6orBgwfjr3/9K9555x389NNPxo6XiMhiFFpZ80It7tJOlszgBobDhg3DsGHDsHXrViQnJ8PZ2Rlnz55FamoqUlNT8d133+Htt99GVVWVMeMlIrIY2g1K/ays4qMd6sotq0Njk9riN18l2yK6c7N2fy4AGDlypG74C0CHw1RERNauULdBqXVVfHxcHOBoL0OdUoWcslr08+ll6pCIOq1b03SJRNKdlyciMmvaio+1NC/UkkgkCG6u+nCeD1ka1ieJiLqJbo6PFS1l1+I8H7JUTHyIiLpJoZVtV9GatokhKz5kaZj4EBF1A7Va0CU+1raqC2DFhywXEx8iom5QUtMIlVqARAJ497K+xIdzfMhSiU585syZg8OHDxsjFiIiq6Fd0eXlbF1dm7W0FZ/s0lo0qdQmjoao80T/bayoqMCECRMwcOBA/POf/0Rubq4x4iIismiF2hVdVtbDR8vfVQG5nRRNagE3ytvftJrI3IhOfPbs2YPc3Fw8//zz+PrrrxESEoIpU6Zgx44dUCqVxoiRiMjiFFhp12YtqVSCYE/NcFcG5/mQBTFK/dXHxwexsbH4/fffcfz4cQwYMACzZs1CYGAgFi9erNfkkIjIFuhWdFlpxQcA+jYnPlmlnOdDlsOoA895eXlITExEYmIiZDIZpk6dirS0NAwZMgRr16415q2IiMyatVd8AKBv8wTnbCY+ZEFEJz5KpRI7d+7EtGnTEBwcjG+//RaLFi3CjRs3sG3bNuzfvx/ffPMNVq5caYx4iYgsgrV2bW5NV/Hhyi6yIKL36goICIBarcbMmTNx4sQJRERE3HLOfffdB3d3d7G3IiKyGEXafbqssGuzlm5JOys+ZEFEJz4LFy7Eyy+/DCcnJ73jgiAgOzsbffv2hbu7OzIyMsTeiojIYuh2ZreBik92aS0EQeD+jGQRRA91vfHGG6iurr7leGlpKUJDQ8VenojI4qjUAoqqrT/x6eOhSXyqG5pQWtNo4miIOkd04iMIQpvHq6uroVBY7194IqL2lOp1bZabOpxuo7CXwb85sePKLrIUBg91xcbGAgAkEglWrFihN9SlUqlw/PjxNuf7EBFZO+2KLi9nB9hZYdfm1vp6OiG/sh5ZpbUY3tfD1OEQ3ZbBic+ZM2cAaCo+aWlpkMtb/qtGLpcjPDwcS5YsER8hEZGF0W5XYc1L2bX6ejnhxPVSruwii2Fw4nPgwAEAwNy5c/Hhhx/C1dXVaEEREVmylu0qbCDxYRNDsjCiV3Vt2bLFGHEQEVkNW1jRpaVNfLiknSyFQYlPbGwsVq1aBWdnZ91cn/bEx8cbFBgRkaUq0PbwsYXEh92bycIYlPicOXNGtwGpdq5PW9jTgYhsUWGl9Tcv1NJWfPIr61GvVEFhLzNxREQdMyjx0c7vufnPRETUskGpLQx1eTnL4SyXoaZRhZyyOgzw7WXqkIg6JHqdZV1dHWprW0qcmZmZ+OCDD7Bv3z6xlyYiski2sEGplkQiQZAnh7vIcohOfGbMmIHt27cDAMrLyxEVFYU1a9ZgxowZWL9+vegAiYgsiUotoKhKu6rL+is+QKsJziU1Jo6E6PZEJz4pKSm45557AAA7duyAv78/MjMzsX37dnz00UeiAyQisiQlNQ1QC7D6rs2taTcrzSqtM3EkRLcnOvGpra2Fi4sLAGDfvn145JFHIJVKMWrUKGRmZooOkIjIkmh7+Hj3sv6uzVrs5UOWRPTfygEDBmDPnj3Izs7Gzz//jIkTJwIACgsL2dSQiGxOgQ2t6NLq6+UMAMgq5VAXmT/Ric+KFSuwZMkShISEYOTIkYiOjgagqf4MHz5cdIBERJbEllZ0abWu+LS3cTWRuRDdufnRRx/FmDFjkJeXh/DwcN3x8ePH4+GHHxZ7eSIii2JLK7q0ers7QioB6pVqFFU12ETjRrJcohMfAPD394e/v7/esaioKGNcmojIomi3q/CxkRVdACC3kyLAzRG55XXIKq1l4kNmzSiJT1JSEpKSklBYWAi1Wq333ubNm41xCyIii1BkQzuzt9bX00mX+Nwd4mnqcIjaJXqOz5tvvomJEyciKSkJxcXFKCsr03sREdkS3QalNlTxAVqWtGeWcGUXmTfRFZ+EhARs3boVs2bNMkY8REQWTbeqy8YqPuzeTJZCdMWnsbERMTExxoiFiMiiqdQCiqttb1UX0Kp7MxMfMnOiE5958+bhyy+/NEYsREQWraRa07VZKtFs3mlLWro3M/Eh8yZ6qKu+vh4bNmzA/v37MWzYMNjb2+u9Hx8fL/YWREQWQTu/x8uGujZraSs+RVUNqGtUwVEuM3FERG0TnficPXsWERERAIBz587pvSeRSMRenojIYhTa6IouAHB3ksNVYYfK+iZkldZikL+LqUMiapPoxOfAgQPGiIOIyOLlVTQnPja2okurr5cTzuVWMvEhs2aUWuyvv/6KJ598EjExMcjNzQUA/Pvf/8Zvv/1mjMsTEVmEnDLN7uR9PBxNHIlp6CY4l3DPLjJfohOfnTt3YtKkSXB0dERKSgoaGjRj3BUVFfjnP/8pOkAiIkuRU6aZ2NvHw8nEkZhGX0/NZqVc0k7mTHTi89ZbbyEhIQEbN27Um9g8evRopKSkiL08EZHFYMWHK7vI/IlOfC5duoR77733luNubm4oLy8Xe3kiIovRkvjYZsVH172ZiQ+ZMdGJj7+/P65evXrL8d9++w39+vUTe3kiIotQr1TpmhfaesUnp7QOarVg4miI2iY68Xn22WexcOFCHD9+HBKJBDdu3MAXX3yBJUuW4PnnnzdGjEREZi+3XFPtcZbL4O5kf5uzrVOAmwJ2UgkaVWrkN2/dQWRuRC9nX7ZsGdRqNcaPH4/a2lrce++9cHBwwJIlS/Diiy8aI0YiIrPXepjLVnuY2cmk6O3hiMySWmSV1iLQ3TYrX2TeRFd8JBIJ/vGPf6C0tBTnzp3DsWPHUFRUhFWrVhkjPiIii9Cyosu2/7HnBGcyd6IrPlpyuRxDhgwx1uWIiCyKtuLTm4kPACCrhIkPmSeDEp/Y2NhOn8u9uojIFtj6UnYtVnzI3BmU+Jw5c0bv65SUFDQ1NWHQoEEAgMuXL0MmkyEyMlJUcIcPH8Z7772H06dPIy8vD7t378ZDDz3U4WcOHjyI2NhYnD9/HkFBQXjttdfw1FNPiYqDiOh2bL15oZauezMTHzJTBiU+rffnio+Ph4uLC7Zt2wYPDw8AQFlZGebOnYt77rlHVHA1NTUIDw/H008/jUceeeS252dkZOCBBx7A/Pnz8cUXXyApKQnz5s1DQEAAJk2aJCoWIqKOsOKj0be5lw+7N5O5Ej3HZ82aNdi3b58u6QEADw8PvPXWW5g4cSJefvllg689ZcoUTJkypdPnJyQkIDQ0FGvWrAEADB48GL/99hvWrl3LxIeIuk29UoWiKm0PH1Z8AKC0phFV9Uq4KGxzaT+ZL9GruiorK1FUVHTL8aKiIlRVVYm9fJckJydjwoQJescmTZqE5OTkdj/T0NCAyspKvRcRUVfcaO7h4ySXwcNGe/houSjs4eksB8B5PmSeRCc+Dz/8MObOnYtdu3YhJycHOTk52LlzJ5555plODU8ZU35+Pvz8/PSO+fn5obKyEnV1dW1+Ji4uDm5ubrpXUFBQT4RKRFak9TCXrfbwaS3Ik8NdZL5EJz4JCQmYMmUK/va3vyE4OBjBwcH429/+hsmTJ+PTTz81Rozdavny5aioqNC9srOzTR0SEVkY3VJ2NuwDAARrJzhzSTuZIdFzfJycnPDpp5/ivffew7Vr1wAA/fv3h7Ozs+jgusrf3x8FBQV6xwoKCuDq6gpHx7Z/ITk4OMDBwaEnwiMiK8UVXfq4pJ3MmdEaGDo7O2PYsGHGupxBoqOj8eOPP+odS0xMRHR0tIkiIiJbwBVd+pj4kDkTPdTVnaqrq5GamorU1FQAmuXqqampyMrKAqAZppo9e7bu/Pnz5yM9PR1Lly7FxYsX8emnn+Kbb77B4sWLTRE+EdkIVnz0aZe0M/Ehc2TWic+pU6cwfPhwDB8+HICmY/Tw4cOxYsUKAEBeXp4uCQKA0NBQ/PDDD0hMTER4eDjWrFmDTZs2cSk7EXUr7c7srPhoaCs+uWV1aFKpTRwNkT6jDXV1h3HjxkEQhHbf37p1a5ufubmzNBFRd2loUqGgUtvDh4kPAPi7KiCXSdGoUiOvol63yovIHJh1xYeIyNzdKK8HADjay3T9a2ydVCpBH09NEsjhLjI3ohOfOXPm4PDhw8aIhYjI4rTM72EPn9b6ckk7mSnRiU9FRQUmTJiAgQMH4p///Cdyc3ONERcRkUXQ9fDhMJeeYK7sIjMlOvHZs2cPcnNz8fzzz+Prr79GSEgIpkyZgh07dkCpVBojRiIis9W64kMt2L2ZzJVR5vj4+PggNjYWv//+O44fP44BAwZg1qxZCAwMxOLFi3HlyhVj3IaIyOy09PDhBN7Wgr00TWwzS2tMHAmRPqNObs7Ly0NiYiISExMhk8kwdepUpKWlYciQIVi7dq0xb0VEZBbYvLBtuiaGnONDZkZ04qNUKrFz505MmzYNwcHB+Pbbb7Fo0SLcuHED27Ztw/79+/HNN99g5cqVxoiXiMis5LLi06ag5lVdlfVNKK9tNHE0RC1E9/EJCAiAWq3GzJkzceLECURERNxyzn333Qd3d3extyIiMisNTSoUVGmWs7Pio89JbgcfFwcUVTUgq7QW7k5c6k/mQXTFZ+HChcjJycG6dev0kh5BEHRdld3d3ZGRkSH2VkREZiWvvB6CACjspfBiD59bcM8uMkeiE5833ngD1dXVtxwvLS1FaGio2MsTEZmt1hOb2cPnVsHs5UNmSHTi096WEtXV1VAoFGIvT0RktrRL2Xu7c5irLVzSTubI4Dk+sbGxAACJRIIVK1bAyallYp9KpcLx48fbnO9DRGQtuKKrY+zeTObI4MRHuxGoIAhIS0uDXN4yvi2XyxEeHo4lS5aIj5CIyEy1NC/kiq62BHtxjg+ZH4MTnwMHDgAA5s6diw8//BCurq5GC4qIyBLklrPi0xFtxSevog6NTWrI7bgvNpme6J/CLVu2MOkhIpvEoa6O+bg4QGEvhVpoSRKJTM2gik9sbCxWrVoFZ2dn3Vyf9sTHxxsUGBGROWtsUiO/UtvDh0NdbZFIJOjr6YTLBdXIKq1FqLezqUMiMizxOXPmjG4DUu1cn7ZweScRWau8ijoIAuBgJ4V3L/bwaY8u8SmpAeBj6nCIDEt8tPN7bv4zEZGt0A5z9fZw5H/kdaCvp6bKwwnOZC4404yIyABc0dU5fZv37GLiQ+bC4Dk+ncU5PkRkjTixuXOCvTQVH/byIXNh8ByfzmD5l4isFROfzmndvVkQBP67QCYneo4PEZEtym21Txe1r4+HIyQSoKZRhZKaRnj3cjB1SGTjOMeHiMgALXN8WPHpiMJehkA3zTNKL6oxcTRE7ONDRNRl+j18mPjczkC/Xsgtr8PlgipEhXqaOhyycezjQ0TURfkV9VA39/Dx4dDNbd3h54KDl4pwpaDK1KEQsY8PEVFXaYe5eruzh09nDPTtBQC4XFBt4kiIjDzHRxAECIJgzEsSEZmd1s0L6fbu8HMBAFwpZMWHTM8oic9nn32GsLAwKBQKKBQKhIWFYdOmTca4NBGR2WHzwq4Z0FzxKa5uRGlNo4mjIVtn0FBXaytWrEB8fDxefPFFREdHAwCSk5OxePFiZGVlYeXKlaKDJCIyJznl7OHTFc4Odujj4YicMs0E51H9vEwdEtkw0YnP+vXrsXHjRsycOVN37MEHH8SwYcPw4osvMvEhIqvD5oVdd4efC3LK6nCFiQ+ZmOihLqVSibvvvvuW45GRkWhqahJ7eSIis8PmhV030I8TnMk8iE58Zs2ahfXr199yfMOGDXjiiSfEXp6IyKwoVWrkVWgSnyBWfDrtDl/NBOdLXNJOJiZ6k1KJRIJNmzZh3759GDVqFADg+PHjyMrKwuzZs40TJRGRmdD28JHbSbn9QhcM8m9e2VVQxT27yKSMsklpZGQkAODatWsAAG9vb3h7e+P8+fMiwyMiMi/ZrXr4SKX8x7uz+vv0gkQClNUqUVzdCB8XJo1kGtyklIioCzix2TCOchn6ejohs6QWVwqqmPiQyYhe1aV14cIFZGVlobGxpUeDRCLB9OnTjXULIiKTY+JjuIG+LsgsqcXlgirEDPA2dThko0QnPunp6Xj44YeRlpYGiUSi69ysHb9VqVRib0FEZDa4ostwd/j1wv4/CnC5kCu7yHREr+pauHAhQkNDUVhYCCcnJ5w/fx6HDx/G3XffjYMHDxohRCIi89HStZkVn67SbV3BlV1kQqIrPsnJyfjll1/g7e0NqVQKqVSKMWPGIC4uDi+99FKHu7cTEVkaDnUZrnUvH67sIlMRXfFRqVRwcdFk8d7e3rhx4wYAIDg4GJcuXRJ7eSIis9GkUiO/sh4Ah7oM0d+nF6QSoKJOiaKqBlOHQzZKdMUnLCwMv//+O0JDQzFy5EisXr0acrkcGzZsQL9+/YwRIxGRWcirqIdKLUAuk8KHPXy6TGEvQ7CXMzKKa3C5oBq+rgpTh0Q2SHTF57XXXoNarQYArFy5EhkZGbjnnnvw448/4qOPPhIdIBGRudAOcwW6K9jDx0ADfbXDXZznQ6YhuuIzadIk3Z8HDBiAixcvorS0FB4eHhy/JSKr0jKxmcNchhrk74J9FwpwpZCJD5mG0fr4ANAtZff09DTmZYmIzEJuOSc2izWweWXXpXwmPmQaooe6AOCzzz5DWFgYFAoFFAoFwsLCsGnTJmNcmojIbHBFl3h3NK/sutK8souop4mu+KxYsQLx8fF48cUXER0dDUCzxH3x4sXIysrCypUrRQdJRGQOONQlXqi3M2RSCaoampBfWY8ANyaR1LNEJz7r16/Hxo0bMXPmTN2xBx98EMOGDcOLL77IxIeIrAYrPuI52MkQ4uWEa0WalV1MfKiniR7qUiqVuPvuu285HhkZiaamJrGXJyIyC00qNfIq2MPHGNjBmUxJdOIza9YsrF+//pbjGzZswBNPPCH28kREZiG/UtPDx14mgS93FhdFO8GZS9rJFAwa6oqNjdX9WSKRYNOmTdi3bx9GjRoFADh+/DiysrIwe/Zs40RJRGRiLT18HNnDR6Q7Wm1dQdTTDEp8bt5/KzIyEgBw7do1AJqtK7y9vXH+/HmR4RERmQfO7zEe7VDX1ULu2UU9z6DE58CBA8aOg4jIrOVqEx93zu8RK8TLGXZSCaobmnCjoh693ZlMUs8xSh+f8vJyrFmzBvPmzcO8efOwdu1aVFRUGOPSWLduHUJCQqBQKDBy5EicOHGi3XO3bt0KiUSi91IouBcMEYnXspSd/0iLJbeTop+PMwDO86GeJzrxOXXqFPr374+1a9eitLQUpaWliI+PR//+/ZGSkiLq2l9//TViY2Px+uuvIyUlBeHh4Zg0aRIKCwvb/Yyrqyvy8vJ0r8zMTFExEBEBrYa6PJn4GMNAruwiExGd+CxevBgPPvggrl+/jl27dmHXrl3IyMjAtGnTsGjRIlHXjo+Px7PPPou5c+diyJAhSEhIgJOTEzZv3tzuZyQSCfz9/XUvPz+/Du/R0NCAyspKvRcR0c1yytm80Jju8NWu7OIEZ+pZRqn4/P3vf4edXct0ITs7OyxduhSnTp0y+LqNjY04ffo0JkyYoDsmlUoxYcIEJCcnt/u56upqBAcHIygoCDNmzLjtBOu4uDi4ubnpXkFBQQbHTETWqUmlRl65tocPKz7G0LKyixUf6lmiEx9XV1dkZWXdcjw7OxsuLi4GX7e4uBgqleqWio2fnx/y8/Pb/MygQYOwefNmfPfdd/j888+hVqsRExODnJycdu+zfPlyVFRU6F7Z2dkGx0xE1qmgqgFNagF2Ugl8XThv0BhahrqqoVZzzy7qOaK3rHj88cfxzDPP4P3330dMTAwA4MiRI3jllVf0trHoCdHR0br9wgAgJiYGgwcPxr/+9S+sWrWqzc84ODjAwYHNyIiofdmlmmGuQHdHyNjDxyhCvJwgl0lRp1Qht7wOQZ4cQqSeITrxef/99yGRSDB79mzdFhX29vZ4/vnn8c477xh8XW9vb8hkMhQUFOgdLygogL+/f6euYW9vj+HDh+Pq1asGx0FEdCarHAAwyN/wKjbps5NpVnZdzK/C5YIqJj7UY0QPdcnlcnz44YcoKytDamoqUlNTUVpairVr14qqpMjlckRGRiIpKUl3TK1WIykpSa+q0xGVSoW0tDQEBAQYHAcRUXJ6CQAgup+XiSOxLi1bV3CCM/UcUYmPUqnE+PHjceXKFTg5OWHo0KEYOnQonJyMk7nHxsZi48aN2LZtG/744w88//zzqKmpwdy5cwEAs2fPxvLly3Xnr1y5Evv27UN6ejpSUlLw5JNPIjMzE/PmzTNKPERke5QqNU5dLwUARPdn4mNMd/hqJjhzSTv1JFFDXfb29jh79qyxYrnF448/jqKiIqxYsQL5+fmIiIjA3r17dROes7KyIJW25G5lZWV49tlnkZ+fDw8PD0RGRuLo0aMYMmRIt8VIRNbtbE45ahtV8HCyxyA/DnUZk67iU8jEh3qORBAEUdPpFy9eDAcHB1HzecxJZWUl3NzcUFFRAVdXV1OHQ0Qm9nHSFaxJvIwpYf5Y/2SkqcOxKhnFNbjv/YNQ2Etx4c3J3PyVROnsv9+iJzc3NTVh8+bN2L9/PyIjI+Hs7Kz3fnx8vNhbEBGZjHZ+TwyHuYyur6cTHOykqFeqkV1Wi2Av59t/iEgk0YnPuXPncNdddwEALl++rPced9wlIkvW0KTC6cwyAJzf0x1kUgn6+/TChbxKXC6oZuJDPUJ04sOd2onIWp3JKkdDkxo+Lg7o79PL1OFYpTv8tIlPFe4f0vEWQ0TGYPCqLrVajXfffRejR4/GiBEjsGzZMtTV1RkzNiIik0q+phnmGtXPixXsbsLNSqmnGZz4vP3223j11VfRq1cv9O7dGx9++CEWLFhgzNiIiExKm/iwf0/3uYO9fKiHGZz4bN++HZ9++il+/vln7NmzB//973/xxRdfQK1WGzM+IiKTqGtU4Uy2Zn4PJzZ3H+1mpVeLqqHinl3UAwxOfLKysjB16lTd1xMmTIBEIsGNGzeMEhgRkSmdziyDUiUgwE2BYC9up9BdgjycoLCXorFJjcySGlOHQzbA4MSnqakJCoX+LsX29vZQKpWigyIiMrXk9GIAmmEuzu/pPlKpBAOaOzhzuIt6gsGrugRBwFNPPaW3H1d9fT3mz5+v18tn165d4iIkIjIB3cRmDnN1uzt8XXAutxJXCqowOaxzm1ATGcrgxGfOnDm3HHvyySdFBUNEZA6qG5rwe04FAE5s7gl3+Gu3rmDFh7qfwYnPli1bjBkHEZHZOHm9FCq1gCBPRwR5cn5Pd9NOcOaSduoJonZnJyKyRse4jL1HDfTVVHzSi2rQpOLKYOpeTHyIiG6i3Z+L21T0jN7ujnCSy9CoUuN6Sa2pwyErx8SHiKiVijolzuVq5/d4mzga2yCVSjDQl8Nd1DOY+BARtXIioxRqAQj1doa/m+L2HyCjGMgOztRDmPgQEbWi26aCw1w9SjvB+XIhKz7UvZj4EBG1opvfw4nNPYqblVJPYeJDRNSsrKYRf+RVAtDsyE49R7tZaXpRDRqaVCaOhqwZEx8iombHMzTVnoG+veDj4nCbs8mYAt0UCHBToEkt4Ke0fFOHQ1aMiQ8RUbOjnN9jMhKJBE+M7AsA2HL0ummDIavGxIeIqJl2YnMMEx+TmBnVF3I7KX7PLseZrDJTh0NWiokPERGAoqoGXCmshkQCjAxl4mMKXr0c8GB4IABgK6s+1E2Y+BARATjWvJrrT/6u8HCWmzga2/VUTAgA4IezeSiorDdtMGSVmPgQEYHL2M1FWG833B3sgSa1gC+OZ5k6HLJCTHyIiMDGhebkqdEhAIAvj2dxaTsZHRMfIrJ5+RX1yCiugVQCRIV6mjocmzfpTn/4uypQXN2AH9PyTB0OWRkmPkRk85LTiwFohlncHO1NHA3Zy6SYFR0MANhy5DoEQTBxRGRNmPgQkc3TDXNxfo/Z+OuIIMjtpDibU4Ez2eWmDoesCBMfIrJ52onNozi/x2x49XLADO3S9iPXTRsMWRUmPkRk07JLa5FdWgeZVIIRIZzfY07mNC9t/zGNS9vJeJj4EJFN01Z7hvVxQy8HOxNHQ62F9XbDiBAubSfjYuJDRDbtGLepMGtPxYQCAL48nsml7WQUTHyIyGYJgtCqcaG3iaOhtky80w8BbgoUVzfih7Nc2k7iMfEhIpt1LrcSeRX1sJdJEBnsYepwqA32MimeHMWl7WQ8THyIyCY1qdR4dXcaAGDinf5wlMtMHBG1R7tre1puBVKyyk0dDlk4Jj5EZJM2/pqBtNwKuCrs8Pq0IaYOhzrg6SxvWdrOXdtJJCY+RGRzrhVVY+3+ywCAFdPvhK+rwsQR0e1ol7b/xKXtJBITHyKyKSq1gKU7zqKxSY2xd/jg/93V29QhUSeE9XZDVIinZmn7sUxTh0MWjIkPEdmU7cnXcTqzDM5yGf75yFBIJBJTh0SdpN21/Qvu2k4iMPEhIpuRXVqL1XsvAQCWTR2M3u6OJo6IumLiEM3S9pKaRnz/O5e2k2GY+BCRTRAEAct2nUWdUoWRoZ54IqqvqUOiLrJrtbR93YGruJhfaeKIyBIx8SEim/D1yWwcuVoChb0U7/6/YZBKOcRliWZG9YW7kz3Si2sw5cNfseTb33GjvM7UYZEFYeJDRFYvr6IOb//wBwBgycRBCPF2NnFEZChPZzm+WzAaDwwNgCAAO07n4L73D+Kdny6iok5p6vDIAjDxISKrJggCXtt9DlUNTYgIcsfc0aGmDolECvZyxron7sLu/y8GUaGeaGhSI+HQNYx97wA2/ZrOic/UISY+RGTVvku9gaSLhZDLpFj96DDIOMRlNYb39cDX/zMKn825GwN9e6G8Vom3fvgDf37/EPacyYVaze0t6FZMfIjIahVVNeCN/54HALz45wG4w8/FxBGRsUkkEowf7IefFt6Dd//fUPi5OiC3vA6Lvk7F9E9+w6HLRdzfi/Qw8SEiq/XG/51Hea0SQwJcMX9cf1OHQ93ITibF4yP64uCS+/DKpEHo5WCH8zcqMWfzCTyWkIwjV4uZABEAJj5EZKX2nsvDD2l5kEklWP3oMNjL+OvOFjjKZVhw3wAcemUcnh4dCrmdFKcyy/DEpuN4fMMxHEsvMXWIZGISgSmwnsrKSri5uaGiogKurq6mDoeIuqisphH/OZmFfx1KR0WdEgvu649XJv3J1GGRieRX1GP9wav4z4lsNKrUAICY/l5YfP8dGBHiaeLoyJg6++83E5+bMPEhskxXC6ux+UgGdqXkoF6p+QcurLcrdsyPgcJeZuLoyNRulNfh04NX8fXJbChVmn/27hnojUUT7kBksIeJoyNjYOJjICY+RJZDEAT8eqUYm49k4OClIt3xIQGueGZMKKaFB8DBjkkPtcgpq8W6A9fw7alsNDWv+hp7hw+eHBWMESEecHeSmzhCMpTVJD7r1q3De++9h/z8fISHh+Pjjz9GVFRUu+d/++23+N///V9cv34dAwcOxLvvvoupU6d2+n5MfIjMX71Shd1ncrH5twxcKawGAEgkwP2D/fD0mFCMDPXk5qPUoezSWnzyy1XsSMmBqtWy90F+LhgR6oERIZ6ICvVEgBv3c7MUVpH4fP3115g9ezYSEhIwcuRIfPDBB/j2229x6dIl+Pr63nL+0aNHce+99yIuLg7Tpk3Dl19+iXfffRcpKSkICwvr1D2Z+BCZF7VaQF5lPdKLqpFeVIMrhVX44Wweymo1XXqd5TL8ZUQQnooJQbAXOzJT12SW1OCz3zLw29VipBfV3PJ+kKejJgkK8cSIUE/0dneEg52UibUZsorEZ+TIkRgxYgQ++eQTAIBarUZQUBBefPFFLFu27JbzH3/8cdTU1OD777/XHRs1ahQiIiKQkJDQqXve/OAamlQoqmowzjdERB0qq1HiWlE10ouqca24BulFNcgortbN2Wmtj4cjnooJwV9GBMFVYW+CaMnaFFc34NT1UpzIKMOJ6yW4cKMSbfVAlEklcJLL4Cy3g5OD5n+dHbRf28FZLoODnRQO9jLIZVLI7aRwsNP8r+bPMs2fZVLI7SSQSaWwl0pgJ5PCTiaBvVTzv3baY1IJBAFQCQLUggBBEKBSA2pBgEottPueWi1ALTSf1+q9rvyzL5VIIJNKIJFovm+ppPnPEgmkzV9Lb35Pd1zznvYagKYyaywSiQS93Vsqcp1NfOyMF4JxNTY24vTp01i+fLnumFQqxYQJE5CcnNzmZ5KTkxEbG6t3bNKkSdizZ0+792loaEBDQ0tiU1mpv9vvpfwqPPjJEQO+AyIyFjupBH29nNDPuxf6+zjjrmAPjP+TL+y4RJ2MyLuXAyaHBWByWAAAoKpeiZSscpzMKMWJ66VIzS5HY5MaKrWAqvomVNU3mThi26awl+Liqild/pzZJj7FxcVQqVTw8/PTO+7n54eLFy+2+Zn8/Pw2z8/Pz2/3PnFxcXjzzTfbfV8CCRzs+MuVqCe4KOzQz7sX+vk4a17Nfw7ydGIfHupxLgp7jL3DB2Pv8AEAqNQCahubUNuoQk1Dy//WNDahpkGF2ub/rWloQqNKjYYmNRqb1GhoUrX6s+Z/tceb1AKUKgFNKjWa1AKa1Go0qZqPNf+5Sa3Wr6BIJZBJJJBIJJBJ0fKetLkSo6vG3PqepPkanSEAUAtorhQJN/25VTVJW3XSHde8p/u6VeXJmAxduGC2iU9PWb58uV6VqLKyEkFBQbqvh/Zxw6W3up5REhGRdZFJJXBR2MOFQ6sWzWwTH29vb8hkMhQUFOgdLygogL+/f5uf8ff379L5AODg4AAHBwfxARMREZHZM9vasVwuR2RkJJKSknTH1Go1kpKSEB0d3eZnoqOj9c4HgMTExHbPJyIiIttithUfAIiNjcWcOXNw9913IyoqCh988AFqamowd+5cAMDs2bPRu3dvxMXFAQAWLlyIsWPHYs2aNXjggQfw1Vdf4dSpU9iwYYMpvw0iIiIyE2ad+Dz++OMoKirCihUrkJ+fj4iICOzdu1c3gTkrKwtSaUvRKiYmBl9++SVee+01vPrqqxg4cCD27NnT6R4+REREZN3Muo+PKbCBIRERkeXp7L/fZjvHh4iIiMjYmPgQERGRzWDiQ0RERDaDiQ8RERHZDCY+REREZDOY+BAREZHNYOJDRERENoOJDxEREdkMJj5ERERkM8x6ywpT0DayrqysNHEkRERE1Fnaf7dvtyEFE5+blJSUAACCgoJMHAkRERF1VVVVFdzc3Np9n4nPTTw9PQFoNkDt6MHRrSorKxEUFITs7Gzuc9ZFfHaG47MzHJ+d4fjsDNddz04QBFRVVSEwMLDD85j43ES727ubmxt/mA3k6urKZ2cgPjvD8dkZjs/OcHx2huuOZ9eZggUnNxMREZHNYOJDRERENoOJz00cHBzw+uuvw8HBwdShWBw+O8Px2RmOz85wfHaG47MznKmfnUS43bovIiIiIivBig8RERHZDCY+REREZDOY+BAREZHNYOJDRERENsMmE59169YhJCQECoUCI0eOxIkTJzo8/9tvv8Wf/vQnKBQKDB06FD/++GMPRWp+uvLstm7dColEovdSKBQ9GK15OHz4MKZPn47AwEBIJBLs2bPntp85ePAg7rrrLjg4OGDAgAHYunVrt8dpjrr67A4ePHjLz5xEIkF+fn7PBGxG4uLiMGLECLi4uMDX1xcPPfQQLl26dNvP8fedYc+Ov+801q9fj2HDhumaE0ZHR+Onn37q8DM9/TNnc4nP119/jdjYWLz++utISUlBeHg4Jk2ahMLCwjbPP3r0KGbOnIlnnnkGZ86cwUMPPYSHHnoI586d6+HITa+rzw7QdObMy8vTvTIzM3swYvNQU1OD8PBwrFu3rlPnZ2Rk4IEHHsB9992H1NRULFq0CPPmzcPPP//czZGan64+O61Lly7p/dz5+vp2U4Tm69ChQ1iwYAGOHTuGxMREKJVKTJw4ETU1Ne1+hr/vNAx5dgB/3wFAnz598M477+D06dM4deoU/vznP2PGjBk4f/58m+eb5GdOsDFRUVHCggULdF+rVCohMDBQiIuLa/P8v/zlL8IDDzygd2zkyJHCc889161xmqOuPrstW7YIbm5uPRSdZQAg7N69u8Nzli5dKtx55516xx5//HFh0qRJ3RiZ+evMsztw4IAAQCgrK+uRmCxJYWGhAEA4dOhQu+fw913bOvPs+PuufR4eHsKmTZvafM8UP3M2VfFpbGzE6dOnMWHCBN0xqVSKCRMmIDk5uc3PJCcn650PAJMmTWr3fGtlyLMDgOrqagQHByMoKKjDrJ9a8GdOvIiICAQEBOD+++/HkSNHTB2OWaioqADQshFzW/iz17bOPDuAv+9uplKp8NVXX6GmpgbR0dFtnmOKnzmbSnyKi4uhUqng5+end9zPz6/dOQD5+fldOt9aGfLsBg0ahM2bN+O7777D559/DrVajZiYGOTk5PREyBarvZ+5yspK1NXVmSgqyxAQEICEhATs3LkTO3fuRFBQEMaNG4eUlBRTh2ZSarUaixYtwujRoxEWFtbuefx9d6vOPjv+vmuRlpaGXr16wcHBAfPnz8fu3bsxZMiQNs81xc8cd2enbhMdHa2X5cfExGDw4MH417/+hVWrVpkwMrJWgwYNwqBBg3Rfx8TE4Nq1a1i7di3+/e9/mzAy01qwYAHOnTuH3377zdShWJzOPjv+vmsxaNAgpKamoqKiAjt27MCcOXNw6NChdpOfnmZTFR9vb2/IZDIUFBToHS8oKIC/v3+bn/H39+/S+dbKkGd3M3t7ewwfPhxXr17tjhCtRns/c66urnB0dDRRVJYrKirKpn/mXnjhBXz//fc4cOAA+vTp0+G5/H2nryvP7ma2/PtOLpdjwIABiIyMRFxcHMLDw/Hhhx+2ea4pfuZsKvGRy+WIjIxEUlKS7pharUZSUlK744/R0dF65wNAYmJiu+dbK0Oe3c1UKhXS0tIQEBDQXWFaBf7MGVdqaqpN/swJgoAXXngBu3fvxi+//ILQ0NDbfoY/exqGPLub8fddC7VajYaGhjbfM8nPXLdNmzZTX331leDg4CBs3bpVuHDhgvA///M/gru7u5Cfny8IgiDMmjVLWLZsme78I0eOCHZ2dsL7778v/PHHH8Lrr78u2NvbC2lpaab6Fkymq8/uzTffFH7++Wfh2rVrwunTp4W//vWvgkKhEM6fP2+qb8EkqqqqhDNnzghnzpwRAAjx8fHCmTNnhMzMTEEQBGHZsmXCrFmzdOenp6cLTk5OwiuvvCL88ccfwrp16wSZTCbs3bvXVN+CyXT12a1du1bYs2ePcOXKFSEtLU1YuHChIJVKhf3795vqWzCZ559/XnBzcxMOHjwo5OXl6V61tbW6c/j7rm2GPDv+vtNYtmyZcOjQISEjI0M4e/assGzZMkEikQj79u0TBME8fuZsLvERBEH4+OOPhb59+wpyuVyIiooSjh07pntv7Nixwpw5c/TO/+abb4Q77rhDkMvlwp133in88MMPPRyx+ejKs1u0aJHuXD8/P2Hq1KlCSkqKCaI2Le0S65tf2mc1Z84cYezYsbd8JiIiQpDL5UK/fv2ELVu29Hjc5qCrz+7dd98V+vfvLygUCsHT01MYN26c8Msvv5gmeBNr67kB0PtZ4u+7thny7Pj7TuPpp58WgoODBblcLvj4+Ajjx4/XJT2CYB4/cxJBEITuqycRERERmQ+bmuNDREREto2JDxEREdkMJj5ERERkM5j4EBERkc1g4kNEREQ2g4kPERER2QwmPkRERGQzmPgQERGRzWDiQ0RERDaDiQ8RERHZDCY+RNRl48aNw6JFi0wdBgDjxzJ48GBs2rSpU+eWlJTA19cX169fN9r9W2v9vZnimf/1r3/FmjVrevSeRN2NiQ+RhXrqqacgkUh0Ly8vL0yePBlnz541dWgWq66uDleuXEF4eHinzn/77bcxY8YMhISEdG9gAHbt2oVVq1Z1+31ae+211/D222+joqKiR+9L1J2Y+BBZsMmTJyMvLw95eXlISkqCnZ0dpk2bZuqwLNa5c+cgCALCwsJue25tbS0+++wzPPPMM+2e09jYaLTYPD094eLiYrTrdUZYWBj69++Pzz//vEfvS9SdmPgQWTAHBwf4+/vD398fERERWLZsGbKzs1FUVAQAaGhowEsvvQRfX18oFAqMGTMGJ0+e1LvGuHHj8NJLL2Hp0qXw9PSEv78/3njjDd37NTU1mD17Nnr16oWAgIBOD33s3bsXY8aMgbu7O7y8vDBt2jRcu3at0/cFgKqqKjzxxBNwdnZGQEAA1q5d2+GQj1qtRlxcHEJDQ+Ho6Ijw8HDs2LHjtrGmpqbiz3/+M8aMGQO1Wo2+ffvigw8+6PAzP/74IxwcHDBq1Ci97+mFF17AokWL4O3tjUmTJnXqWQC3f843f9+duWZnnvGOHTswdOhQODo6wsvLCxMmTEBNTY3u/enTp+Orr7663SMkshhMfIisRHV1NT7//HMMGDAAXl5eAIClS5di586d2LZtG1JSUjBgwABMmjQJpaWlep/dtm0bnJ2dcfz4caxevRorV65EYmIiAOCVV17BoUOH8N1332Hfvn04ePAgUlJSbhtPTU0NYmNjcerUKSQlJUEqleLhhx+GWq3u1H0BIDY2FkeOHMH//d//ITExEb/++muH946Li8P27duRkJCA8+fPY/HixXjyySdx6NChdj9z7do1jB07Fn/+85/x4IMP4pFHHsHLL7+MxYsXIzU1td3P/frrr4iMjLzl+LZt2yCXy3HkyBEkJCR0+ll09Tl35praeNp7xnl5eZg5cyaefvpp/PHHHzh48CAeeeQRCIKg+3xUVBROnDiBhoaGdmMhsigCEVmkOXPmCDKZTHB2dhacnZ0FAEJAQIBw+vRpQRAEobq6WrC3txe++OIL3WcaGxuFwMBAYfXq1bpjY8eOFcaMGaN37REjRgh///vfhaqqKkEulwvffPON7r2SkhLB0dFRWLhwYZfiLSoqEgAIaWlpt72vIAhCZWWlYG9vL3z77be698vLywUnJye9e48dO1ZYuHChUF9fLzg5OQlHjx7Vu+YzzzwjzJw5s924JkyYIDz11FOCIAhCVFSUsGbNGkGlUgmurq7CRx991O7nZsyYITz99NN6x8aOHSsMHz68g6egcfOz6Mxz1n6fnb2m9jMdPePTp08LAITr16+3e93ff//9tucQWRJWfIgs2H333YfU1FSkpqbixIkTmDRpEqZMmYLMzExcu3YNSqUSo0eP1p1vb2+PqKgo/PHHH3rXGTZsmN7XAQEBKCwsxLVr19DY2IiRI0fq3vP09MSgQYNuG9uVK1cwc+ZM9OvXD66urroJwFlZWbe9LwCkp6dDqVQiKipK976bm1u797569Spqa2tx//33o1evXrrX9u3bbxkC0srPz8cvv/yC+fPnQ6VSIS0tDREREZBKpZDJZJDL5e1+f3V1dVAoFLccb6sKdLtnYchz7szzBTp+xuHh4Rg/fjyGDh2Kxx57DBs3bkRZWZne+Y6OjgA0c5qIrIGdqQMgIsM5OztjwIABuq83bdoENzc3bNy4EX/5y186fR17e3u9ryUSyS1DJl01ffp0BAcHY+PGjQgMDIRarUZYWJjehF9j3re6uhoA8MMPP6B379567zk4OLT5mWPHjkGtViMiIgKXLl1CXV0dIiIicP36dZSVlSEmJqbd+3l7e9+SJACa/09u1pln0VWdvWZHz1gmkyExMRFHjx7Fvn378PHHH+Mf//gHjh8/jtDQUADQDYv6+PgYHCuROWHFh8iKSCQSSKVS1NXVoX///rq5JlpKpRInT57EkCFDOnW9/v37w97eHsePH9cdKysrw+XLlzv8XElJCS5duoTXXnsN48ePx+DBg9tMEjrSr18/2Nvb603GrqioaPfeQ4YMgYODA7KysjBgwAC9V1BQUJuf0SYJ9fX1OHPmDIKDg+Hp6YmEhASEhYVh6NCh7cY3fPhwXLhw4bbfR2eeRVefszGer5ZEIsHo0aPx5ptv4syZM5DL5di9e7fu/XPnzqFPnz7w9vY26PpE5oYVHyIL1tDQgPz8fACafyg/+eQTVFdXY/r06XB2dsbzzz+PV155BZ6enujbty9Wr16N2traDpdgt9arVy8888wzeOWVV+Dl5QVfX1/84x//gFTa8X8zeXh4wMvLCxs2bEBAQACysrKwbNmyLn1vLi4umDNnji5+X19fvP7665BKpZBIJG2ev2TJEixevBhqtRpjxoxBRUUFjhw5AldXV8yZM+eWz0RHR8POzg4rV65EdXU1+vXrh08++QQff/wxDh8+3GF8kyZNwvLly1FWVgYPDw9Rz6Krz9kYzxcAjh8/jqSkJEycOBG+vr44fvw4ioqKMHjwYN05v/76KyZOnNjlaxOZKyY+RBZs7969CAgIAKD5h/9Pf/oTvv32W4wbNw4A8M4770CtVmPWrFmoqqrC3XffjZ9//rnDf6hv9t577+mSKRcXF7z88su3bWgnlUrx1Vdf4aWXXkJYWBgGDRqEjz76SBdXZ8XHx2P+/PmYNm0aXF1dsXTpUmRnZ7c5twYAVq1aBR8fH8TFxSE9PR3u7u6466678Oqrr7Z5flBQEDZv3oy///3vyMvLg52dHWpra7F379425+q0NnToUNx111345ptv8Nxzz7V7XmefRVees7Ger6urKw4fPowPPvgAlZWVCA4Oxpo1azBlyhQAmkrYnj17sHfv3i5dl8icSQSh1bpFIiIzVlNTg969e2PNmjWdrlp1lqenJ7Zu3YoHH3yw05/54Ycf8Morr+DcuXO3rYJZovXr12P37t3Yt2+fqUMhMhpWfIjIbJ05cwYXL15EVFQUKioqsHLlSgDAjBkzjHqfnJwclJWVdapjc2sPPPAArly5gtzc3HbnEVkye3t7fPzxx6YOg8ioWPEhIrN15swZzJs3D5cuXYJcLkdkZCTi4+M7nHRsiJ9++gmPPfYYqqqq2pw/RETWg4kPERER2QzrG5QmIiIiagcTHyIiIrIZTHyIiIjIZjDxISIiIpvBxIeIiIhsBhMfIiIishlMfIiIiMhmMPEhIiIim8HEh4iIiGwGEx8iIiKyGf8/8KupvEzErYAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from ase.io import read\n",
    "from matscipy.angle_distribution import angle_distribution\n",
    "from matscipy.neighbours import neighbour_list\n",
    "\n",
    "# Load a disordered configuration (amorphous carbon)\n",
    "a = read('../../tests/aC.cfg')\n",
    "\n",
    "# Two atoms are bonded if their distance is smaller than this value\n",
    "max_bond_distance = 2.0  # Angstroms\n",
    "\n",
    "# Get distance between atoms up to a distance of 2 Angstroms which includes only nearest neighbors.\n",
    "# The capital 'D' returns distance vectors, which are needed for angles.\n",
    "i, j, d = neighbour_list('ijD', a, cutoff=max_bond_distance)\n",
    "\n",
    "# Get compute angle distribution\n",
    "nbins = 50\n",
    "bin_spacing = np.pi/nbins\n",
    "angle = (np.arange(nbins) + 0.5) * bin_spacing\n",
    "nb_angles = angle_distribution(i, j, d, nbins)  # `angle_distribution` returns the number of triangles\n",
    "probability = nb_angles / (np.sum(nb_angles) * bin_spacing)  # Normalize\n",
    "\n",
    "# Plot the angle distribution\n",
    "plt.plot(angle, probability)\n",
    "plt.xlabel(r'Bond angle $\\phi$ (radians)')\n",
    "plt.ylabel(r'Probability density $f_\\text{a}(\\phi)$ (radians$^{-1}$)')\n",
    "plt.xlim(0, np.pi)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.1"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
